package Core::Ascore;
use utf8;
use strict;
use warnings;

our $VERSION = '1.00';
our @EXPORT_OK = qw(ascore convert_sty);
use base 'Exporter';

use Math::Combinatorics;
use Math::NumberCruncher;
use Data::Dumper;

BEGIN {push @INC, "../"}

use Core::IntConsts qw(get_masas_aa get_masas_mod);
use Parsers::Omssa;
use Parsers::Sequest;
use Parsers::Phenyx;

my %masas_aa = %{&get_masas_aa};
my %masas_mod = %{&get_masas_mod};


sub ascore{
	my ($peptide, $parent_mass, $z, $masses,
		$intensities, $precision, $limit_mz_sup, $limit_mz_inf) = @_;
	#El numero de sitios de fosforilacion ha de ser mayor que el de
	#fosforilaciones para que el Ascore tenga sentido.
	#Si el Ascore no tiene sentido, lo devuelvo como 1000
	my $num_fosfo_sites = 0;
	my $num_deshidrat_sites = 0;
	my $num_modif_sites = 0;
	my $num_fosfos = 0;
	my $num_deshidrats = 0;
	my $num_modifs = 0;

	my $pept_conteo_modifs = $peptide;
	while($pept_conteo_modifs =~ s/(\(21\)|\(21\,|21\))/x/){
		$num_fosfos++;
	}
	while($pept_conteo_modifs =~ s/(\(23\)|\(23\,|23\))/x/){
		$num_deshidrats++;
	}
	while($pept_conteo_modifs =~ s/^Y/x/){
		$num_fosfo_sites++;
		$num_modif_sites++;
	}
	while($pept_conteo_modifs =~ s/^[ST]/x/){
		$num_fosfo_sites++;
		$num_deshidrat_sites++;
		$num_modif_sites++;
	}
	while($pept_conteo_modifs =~ s/Y(\w)/x$1/){
		$num_fosfo_sites++;
		$num_modif_sites++;
	}
	while($pept_conteo_modifs =~ s/[ST](\w)/x$1/){
		$num_fosfo_sites++;
		$num_deshidrat_sites++;
		$num_modif_sites++;
	}

	my %result;
	#Si el numero de sitios de modificacion = numero de modificaciones
	$num_modifs = $num_fosfos + $num_deshidrats;
	if($num_modif_sites == $num_modifs){
		$result{'phosphorylations'} = $num_fosfos;
		$result{'dehydratations'} = $num_deshidrats;
		$result{'peptide'} = $peptide;
		if($num_modifs == 1){
			$result{'ascore'} = [1000];
		}
		elsif($num_modifs == 2){
			$result{'ascore'} = [1000, 999];
		}
		elsif($num_modifs == 3){
			$result{'ascore'} = [1000, 999, 998];
		}
		return %result;
	}
	#Caso especial: numero de deshidrataciones = numero posibles sitios de deshidratacion
	if($num_deshidrats > 0){
		if($num_deshidrat_sites == $num_deshidrats){
			if($num_modifs == 1){
				$result{'phosphorylations'} = $num_fosfos;
				$result{'dehydratations'} = $num_deshidrats;
				$result{'peptide'} = $peptide;
				$result{'ascore'} = [1000];
				return %result;
			}
			else {
				#Caso no resuelto: num.DHs == num.sitios.DH Y num_mods>=2
				$result{'phosphorylations'} = $num_fosfos;
				$result{'dehydratations'} = $num_deshidrats;
				$result{'peptide'} = $peptide;
				if($num_modifs == 2){
					$result{'ascore'} = [0, 0];
				}
				elsif($num_modifs == 3){
					$result{'ascore'} = [0, 0, 0];
				}
				print "CASO NO RESUELTO \n";
				return %result;
			}
		}
	}
	if($num_modifs > 0){
		%result = process_ascore($peptide, $precision, $masses,
								 $intensities, $z, $limit_mz_sup, $limit_mz_inf);
	}
	else{
		$result{'peptide'} = $peptide;
	}
	$result{'phosphorylations'} = $num_fosfos;
	$result{'dehydratations'} = $num_deshidrats;
	return %result;
}


sub process_ascore{
	my $peptido = $_[0];        #en formato [styuv]
	my $precision = $_[1];      #en unidades m/z
	my $masses = $_[2];         #referencia a array de masas
	my $intensities = $_[3];    #referencia a array intensidades
	my $z = $_[4];              #carga
	my $limit_mz_sup = $_[5];
	my $limit_mz_inf = $_[6];
	my @masas = split /\|/, $masses;
	my @intensidades = split /\|/, $intensities;
	my @pepts_posibles = peptide_collection($peptido);

	## VARIABLES USADAS
	# @pepts_by_score <- secuencias en formato original ordenadas por score de + a -
	# $best_pept <-----  secuencia peptido de mayor ascore
	# @complements  <--  $site;
	# @col_peptidos <--  [@peptide_aas];   AoA de "peptido trozeados" formato stuv
	# @col_psites   <--  [@temp];     	   AoA de sitios de fosforilacion (incluido pept referencia)
	# @psite_pairs  <--  parejas compltarias para cada psite guardadas como indices de col_psites;
	# @best_psites                         sitios de fosfo en pept de mayor score
	# @levels_max_diff	<---  levels (10) de maxima diferencia de scores entre el primero y todos
	# $num_psites = $#best_psites+1;
	# $num_secuencias = $#col_psites+1;

	#Obtencion de los Peptide Scores
	my %pept_scores;
	my %matches_ascore;
	foreach my $pept(@pepts_posibles){
		my %fragmentos_teoricos = %{prever_fragmentos($pept, $z, $limit_mz_sup, $limit_mz_inf)};
		$matches_ascore{$pept} = matches_ascore(\%fragmentos_teoricos, \@masas,
											   \@intensidades, $precision, $pept);
	}

	my %weigths = (1=>0.5, 2=>0.75, 3=>1, 4=>1, 5=>1, 6=>1, 7=>0.75, 8=>0.5, 9=>0.25, 10=>0.25);
	my $suma_niveles = 7;     #suma de 0.5+0.75+1+1+1+1+0.75+0.5+0.25+0.25

	foreach my $pept(keys %matches_ascore){
		my $suma_datos_matches = 0;
		my $score = 0;
		for (my $lvel=1; $lvel<11; $lvel++){
			$score += $weigths{$lvel} * $matches_ascore{$pept}{$lvel}{score};
			$matches_ascore{$pept}{$lvel}{score_ponderado} =
				($weigths{$lvel} * $matches_ascore{$pept}{$lvel}{score} / $suma_niveles);
		}
		$pept_scores{$pept} = $score / $suma_niveles;
	}

	#Ordeno peptidos segun Peptide score (de mayor a menor)
	my @pepts_by_score = sort {$pept_scores{$b} <=> $pept_scores{$a}} keys %pept_scores;
	my $best_pept = $pepts_by_score[0];

	#Convierto peptidos en formato [styuv]
	my @col_peptidos;
	my @peptide_aas;
	foreach (@pepts_by_score){
		@peptide_aas = split //, convert_sty($_);
		push @col_peptidos, [@peptide_aas];
	}

	#Detecto puntos de modificacion
	my @col_psites;
	foreach my $pept (@col_peptidos){
		my @temp;
		for(0..$#{$pept}){
			if(@{$pept}[$_] =~ /[styuv]/){
				push @temp, $_;
			}
		}
		push @col_psites, [@temp];
	}

	#Busco Peptidos-2 para cada modificacion: Se busca un peptido con el mayor
	# score posible y que mantiene todos los puntos de fosforilacion del
	# peptido de referencia excepto la fosforilacion estudiada.
	#De otra forma, se buscan peptidos con ese punto cambiado
	# pero que el numero de coincidencias sea = total_psites-1
	my @best_psites = @{$col_psites[0]};
	my $num_psites = $#best_psites + 1;
	my $num_secuencias = $#col_psites + 1;
	my @psite_pairs;

	my $psites_of_i;                                    # puntos de fosforilacion de la secuencia i
	for (my $site=0; $site<$num_psites; $site++){       #para cada modificacion en el peptido de referencia
		for (1..$num_secuencias){                       #para cada una de las posibles secuencias
			$psites_of_i = $col_psites[$_];             #obten sus psites
			#mira si el sitio de referencia no esta fosforilado en la posible pareja
			if (!($best_psites[$site] ~~ @{$psites_of_i})){
				my $count = 0;                          #si es asi,
				foreach my $my_psite (@best_psites){    #calcula cuantos coinciden...
					if($my_psite ~~ @{$psites_of_i}){
						$count += 1;
					}
				}
				if ($count == $num_psites - 1){         #si coinciden todos menos 1
					push @psite_pairs, $_;              #guarda el idx de la secuencia pareja
					last;                               #y pasa a otro punto de fosforilacion
				}
			}
		}
	}

	#Obtengo niveles maxima diferencia entre Peptide Scores
	#(P1 respecto a peptido2_f1, peptido2_f2 y peptido2_f3)
	my @levels_max_diff;
	my %datos_bp = %{$matches_ascore{$best_pept}};
	foreach my $pept (@pepts_by_score){
		my $diferencia = 0;
		my $level;
		for (1..10) {
			my $dif = $datos_bp{$_}{score} - ${$matches_ascore{$pept}}{$_}{score};
			if($dif > $diferencia){
				$diferencia = $dif;
				$level = $_;
			}
		}
		if ($diferencia == 0){         #Cuando el pep1 y el pep2 son identicos;
			$level = 5;                # uso el 5 de forma arbitraria
		}
		push @levels_max_diff, $level;
	}

	#Encuentro puntos complementarios a las f1,f2 y f3 en cada peptido
	#Punto complementario es el aminoacido que este modificado en p2 pero NO en p1
	my @complements;
	for(my $i=0; $i < $num_psites; $i++){                      #de psite en psite del peptido referencia
		my @couple_sites = @{$col_psites[$psite_pairs[$i]]};   #cojo los sites del peptido pareja
		foreach my $site (@couple_sites){                      #busco un site en el peptido pareja correspondiente
			if (!($site ~~ @best_psites) && ($site != $best_psites[$i])){   #que no este en la referencia
				#print "para site $i el complementario es $site\n";
				push @complements, $site;
				last;
			}
		}
	}

	#Calculo fragmentos
	my $pept_compl;
	my %fingerp_ref =  %{prever_fragmentos($pepts_by_score[0], $z,
										   $limit_mz_sup, $limit_mz_inf)};
	my @frag_ref;
	my @frag_compl;
	my %fingerp_compl;
	my ($start, $end);
	my $idx = 0;
	foreach my $compsite (@complements){

		if ($best_psites[$idx] > $compsite){
			$start = $compsite + 1;
			$end = $best_psites[$idx] + 1;
		}
		else {
			$start = $best_psites[$idx] + 1;
			$end = $compsite + 1;
		}

		$pept_compl = $pepts_by_score[$psite_pairs[$idx]];

		#print "peptido referencia ---> $pepts_by_score[0]\n";
		#print "peptido_compl $idx ---> $pept_compl\n\n";
		#print "start-end, $start, $end\n";

		%fingerp_compl =  %{prever_fragmentos($pept_compl, $z, $limit_mz_sup, $limit_mz_inf)};

		push @frag_ref, get_specific_frags(convert_sty($pepts_by_score[0]),
												\%fingerp_ref, $start, $end);

		push @frag_compl, get_specific_frags(convert_sty($pept_compl),
												\%fingerp_compl, $start, $end);
		$idx++;
	}

	#Obtengo datos-ascore: con matches de los fragmentos site-determining
	my @ascores;
	for (my $psite=0; $psite<$num_psites; $psite++){
		$pept_compl = $pepts_by_score[$psite_pairs[$psite]];
		my $level = $levels_max_diff[$psite_pairs[$psite]];
		#
		my %dat_ascore_bp = %{matches_ascore($frag_ref[$psite], \@masas,
											  \@intensidades, $precision, $best_pept)};
		my %dat_ascore_compl = %{matches_ascore($frag_compl[$psite], \@masas,
												\@intensidades, $precision, $pept_compl)};
		#
		my $value = $dat_ascore_bp{$level}{score} -
					$dat_ascore_compl{$level}{score};

		if ($value < 0.1){$value = 0;}
		push @ascores, $value;
	}

	my %result;
	$result{'peptide'} = $best_pept;
	$result{'peptide_score'} = $pept_scores{$best_pept};
	$result{'ascore'} = [@ascores];

	return %result;
}


sub peptide_collection{
	my $pept = shift;
	#Limitado a 1/dos chemtags por aa o 1 chemtag + 1 phos
	#			2/Los chemtag no cambian de sitio
	my @posibles_fosfo_sites;
	my @posibles_deshidrat_sites;
	my $num_fosfos = 0;
	my $num_deshidrats = 0;
	my $num_chemtag = 0;
	my $alrdy_one;
	my (@aas, @mods);
	my $aa_num = 0;
	my $aa;
	while ($pept){
		if ($pept =~ s/^(\w)\((\w+),(\w+)\)//){    #Limitado a 2 chemtags por aa.
			$alrdy_one = 0;
			$aa = $1;
			my @mod1 = ($2, $3);
			push @aas, $aa;
			#Calculo el numero de fosfos/chemtags/DH
			foreach (@mod1) {
				if ($_ == 21) {
					$num_fosfos++;
				}
				elsif ($_ == 23) {
					$num_deshidrats++;
				}
				else {
					if ($alrdy_one){
						my $frst = pop @mods;
						push @mods, "$frst,$_";
					} else {
						push @mods, $_;
						$alrdy_one = 1;
					}
					if($aa_num > 0){        #Si hay chemtag en un aa no N-terminal,
						$aa_num++;          #no se puede modificar con nada mas
						next;               #ojo este next tiene que ser un label
					}
				}
			}
		}
		elsif ($pept =~ s/^(\w)\((\w+)\)//){
			$aa = $1;
			push @aas, $aa;
			#Calculo el numero de fosfos/chemtags/DH
			if ($2 == 21){
				push @mods, '';
				$num_fosfos++;
			}
			elsif ($2 == 23){
				push @mods, '';
				$num_deshidrats++;
			}
			else {
				push @mods, $2;
				$num_chemtag++;
				if ($aa_num > 0){        #Si hay chemtag en un aa no N-terminal,
					$aa_num++;           #no se puede modificar con nada mas
					next;
				}
			}
		}
		else {
			$pept =~ s/^(\w)//;
			$aa = $1;
			push @aas, $aa;
			push @mods, '';
		}
		if($aa =~ /[S,T,Y]/){
			push @posibles_fosfo_sites, $aa_num;
			if($aa =~ /[S,T]/){
				push @posibles_deshidrat_sites, $aa_num;
			}
		}
		$aa_num++;
	}
	#Obtengo las diferentes combinaciones de fosforilados y deshidratados
	my @combins_fosfo;
	my @combins_deshidrats;

	my $combin_1 = Math::Combinatorics->new(count => $num_fosfos,
											 data => [@posibles_fosfo_sites]);
	while(my @combo1 = $combin_1->next_combination){
		push @combins_fosfo, \@combo1;
	}
	my $combin_2 = Math::Combinatorics->new(count => $num_deshidrats,
											 data => [@posibles_deshidrat_sites]);

	while(my @combo2 = $combin_2->next_combination){
		push @combins_deshidrats, \@combo2;
	}

	my @combins_peptidos;
	if(@combins_fosfo and @combins_deshidrats){
		#print "Hay phospho y DH\t";
		foreach my $fos(@combins_fosfo){
			my %usados;
			my @fos = @$fos;
			foreach my $f(@fos){
				$usados{$f} = 1;
			}
			LABEL1:
			foreach my $deshidrat(@combins_deshidrats){
				my @deshidrats = @$deshidrat;
				foreach my $d(@deshidrats){
					if(defined $usados{$d}){
						next LABEL1;
					}
				}
				# Recupero chemtags y otras PTMs
				my @mods_new;
				for (0..$#mods){
					if($mods[$_]){
						$mods_new[$_] = "($mods[$_])";
					}
				}
				foreach (@fos){
					if($mods_new[$_]){
						my $t = $mods_new[$_];
						$t =~ s/[\(,\)]//g;
						$mods_new[$_] = "($t,21)";
					}
					else{
						$mods_new[$_]= '(21)';
					}
				}
				foreach (@deshidrats){
					if($mods_new[$_]){
						my $t = $mods_new[$_];
						$t =~ s/[\(,\)]//g;
						$mods_new[$_] = "($t,23)";
					}
					else{
						$mods_new[$_] = '(23)';
					}
				}
				#Genero el peptido
				my $pept_string = get_pept_string(\@aas, \@mods_new);
				push @combins_peptidos, $pept_string;
			}
		}
	}
	elsif (@combins_fosfo){
		#print "solo fosfo\t";
		foreach my $fos(@combins_fosfo){
			my @fos = @$fos;
			# Recupero chemtags y otras PTMs
			my @mods_new;
			for(0..$#mods){
				if($mods[$_]){
					$mods_new[$_] = "($mods[$_])";
				}
			}
			foreach (@fos){
				if($mods_new[$_]){
					my $t = $mods_new[$_];
					$t =~ s/[\(,\)]//g;
					$mods_new[$_] = "($t,21)";
				}
				else{
					$mods_new[$_] = '(21)';
				}
			}
			#Genero el peptido
			my $pept_string = get_pept_string(\@aas, \@mods_new);
			push @combins_peptidos, $pept_string;
		}
	}
	elsif (@combins_deshidrats){
		#print "solo DH\t";
		foreach my $deshidrat(@combins_deshidrats){
			my @deshidrats = @$deshidrat;
			# Recupero chemtags y otras PTMs
			my @mods_new;
			for(0..$#mods){
				if($mods[$_]){
					$mods_new[$_] = "($mods[$_])";
				}
			}
			foreach (@deshidrats){
				if($mods_new[$_]){
					my $t = $mods_new[$_];
					$t =~ s/[\(,\)]//g;
					$mods_new[$_] = "($t,23)";
				}
				else{
					$mods_new[$_] = '(23)';
				}
			}
			#Genero el peptido
			my $pept_string = get_pept_string(\@aas, \@mods_new);
			push @combins_peptidos, $pept_string;
		}
	}
	else{
		push @combins_peptidos, $_[0];
	}
	return @combins_peptidos;
}

sub get_pept_string{
	my @aas = @{$_[0]};
	my @mods_new = @{$_[1]};
	my $pept_string = '';
	for(0..$#aas){
		if($mods_new[$_]){
			$pept_string = $pept_string . $aas[$_] . $mods_new[$_];
		}
		else{
			$pept_string = $pept_string . $aas[$_];
		}
	}
	return $pept_string
}                                                                                    ####226 lineas


sub prever_fragmentos{
	#Predict the spectrum of a sequence. Allows phosphorilations on
	# STY (+80:s,t,y), water loss on ST (-18:u,v), and Met oxidation (+16:m).
	#Calculates y and b series for a given charge state,
	#ARGUMENTS:
	#     - $_[0]-> sequence
	#     - $_[1]-> charge state (1,2 or 3)
	#     - $_[2]-> mz upper limit for fragment generation
	#RETURN: hash{label}->mass.
	my $peptide = $_[0];
	my $z = $_[1];
	my $limit_mz_sup = $_[2];
	my $limit_mz_inf = $_[3];
	my @masas_aas;

	while($peptide){
		if($peptide =~ s/^(\w)\((\w+),(\w+)\)//){
			push @masas_aas, ($masas_aa{$1} + $masas_mod{$2} + $masas_mod{$3});
		}
		elsif($peptide =~ s/^(\w)\((\w+)\)//){
			push @masas_aas, ($masas_aa{$1} + $masas_mod{$2});
		}
		else{
			$peptide =~ s/^(\w)//;
			push @masas_aas, ($masas_aa{$1});
		}
	}
	my $H = 1.007825;
	my $M = 19.018465;              #19.018 = H20 + H
	my $H2O = 18.010565;
	my $aan = $#masas_aas;
	my (@serie_y_1, @serie_b_1, @serie_y_2, @serie_b_2, @serie_y_3, @serie_b_3);

	for (0..$aan){
		$serie_b_1[$_] = $M + $masas_aas[$_] - $H2O;
		$M += $masas_aas[$_];
		}

	for (0..$aan){
		$serie_y_1[$_] = $M - $serie_b_1[($aan - $_)] + $H;
		}

	#La carga maxima permitida es 1 para z=1 o 2,
	if($z > 2){
		for (0..$aan){
			$serie_y_2[$_] = ($serie_y_1[$_] + $H) / 2;
			$serie_b_2[$_] = ($serie_b_1[$_] + $H) / 2;
			}
		if($z > 3){
			for (0..$aan){
				$serie_y_3[$_] = ($serie_y_1[$_] + $H) / 3;
				$serie_b_3[$_] = ($serie_b_1[$_] + $H) / 3;
			}
		}
	}

	my %fragmentos;
	for (0..$aan){
		my $j = $_ + 1;
		if ($_ < ($aan)){                             # El b(n) bo existe
			my $num = $j;
			if($num < 10){$num = '0' . $num};
			$fragmentos{"b$num+1"} = sprintf("%.3f", $serie_b_1[$_]);
			if ($z > 2){
				$fragmentos{"b$num+2"} = sprintf("%.3f", $serie_b_2[$_]);
				if ($z > 3){
					$fragmentos{"b$num+3"} = sprintf("%.3f", $serie_b_3[$_]);
				}
			}
		}
		if(($_ > 0)){										# El y(0) bo existe
			my $num = $_;
			if($num < 10){$num = '0' . $num}
			$fragmentos{"y$num+1"} = sprintf("%.3f", $serie_y_1[$_]);
			if ($z > 2){
				$fragmentos{"y$num+2"} = sprintf("%.3f", $serie_y_2[$_]);
				if($z > 3){
					$fragmentos{"y$num+3"} = sprintf("%.3f", $serie_y_3[$_]);
				}
			}
		}
	}
	#Solo se consideran los fragmentos entre limit_mz_sup e _inf
	foreach (keys %fragmentos){
		if(($fragmentos{$_} > $limit_mz_sup) or ($fragmentos{$_} < $limit_mz_inf)){
			delete $fragmentos{$_};
		}
	}
	return \%fragmentos;
}


sub binomial{
	# ($n,$k,$p) -> $k successes in $n tries, given a probability of $p
	my ($n, $k, $p) = @_;
	my $binomial_acumulada = 0;
	for ($k..$n){
		$binomial_acumulada += Math::NumberCruncher::Binomial($n, $_, $p);
	}
	return $binomial_acumulada;
}


sub convert_sty{
	my $pept = shift;
	#cambio fosfo o DH simples
	$pept =~ s/S\(21\)/s/g;
	$pept =~ s/T\(21\)/t/g;
	$pept =~ s/Y\(21\)/y/g;
	$pept =~ s/S\(23\)/u/g;
	$pept =~ s/T\(23\)/v/g;
	#cambio fosfo o DH combinados con chemtag en Nt
	$pept =~ s/S\((214|737)\,21\)/s/;
	$pept =~ s/T\((214|737)\,21\)/t/;
	$pept =~ s/Y\((214|737)\,21\)/y/;
	$pept =~ s/S\(21\,(214|737)\)/s/;
	$pept =~ s/T\(21\,(214|737)\)/t/;
	$pept =~ s/Y\(21\,(214|737)\)/y/;
	$pept =~ s/S\((214|737)\,23\)/u/;
	$pept =~ s/T\((214|737)\,23\)/v/;
	$pept =~ s/S\(23\,(214|737)\)/u/;
	$pept =~ s/T\(23\,(214|737)\)/v/;
	#Elimino todas las demas modificaciones
	$pept =~ s/\(\d+\)//g;
	return $pept;
}


sub get_specific_frags{
	#Devuelvo los fragmentos que estan
	#entre el primer y ultimo aa que se pueda fosforilar.
	my $peptido = $_[0];
	my %fragmentos = %{$_[1]};
	my $first_aa_fosfable = $_[2];
	my $last_aa_fosfable = $_[3];
	my %fragmentos_filtrados;
	my $len_pept = length $peptido;

	foreach my $frag (keys %fragmentos){
		my $frag2 = (split /\+/, $frag)[0];
		if ($frag2 =~ s/y//){
			#para las y, valen los fragmentos entre [largo-ultimo+1 , largo-primero]
			if(($frag2 >= ($len_pept - $last_aa_fosfable + 1)) and
			   ($frag2 <= ($len_pept - $first_aa_fosfable))){
				$fragmentos_filtrados{$frag} = $fragmentos{$frag};
			}
		}
		elsif ($frag2 =~ s/b//){
			#para las b, valen los fragmentos entre [primero , ultimo-1]
			if(($frag2 >= $first_aa_fosfable) and ($frag2 <= ($last_aa_fosfable-1))){
				$fragmentos_filtrados{$frag} = $fragmentos{$frag};
			}
		}
	}
	%fragmentos = ();
	%fragmentos = %fragmentos_filtrados;    #Uso solo fragmentos en zona de interes
	return \%fragmentos;
}


sub matches_ascore{
	my %fragmentos = %{$_[0]};
	my @masas= @{$_[1]};
	my @intensidades = @{$_[2]};
	my $fragment_error = $_[3];
	my $peptido = $_[4];

	my $num_fragments = keys %fragmentos;
	my $picos_x_sector;
	my %datos_num_picos;
	my @rangos_masas;

	my ($rango_masas_inf, $rango_masas_sup) = (200, 2000);
	for(my $i=$rango_masas_inf; $i<$rango_masas_sup; $i=$i+100){
		push @rangos_masas, $i;
	}

	for ($picos_x_sector=1; $picos_x_sector<11; $picos_x_sector++){
		my @masas_exp_match_intervalo;
		my @masas_teor_match_intervalo;
		my @tags_match_intervalo;

		my $correspondencias = 0;
		my $num_peaks_intervalo = 0;

		for my $intervalo(@rangos_masas){
			my %picos_intervalo;
			my $intervalo_inferior = $intervalo - 100;
			for (0..$#masas){
				if(($masas[$_] > $intervalo_inferior) and ($masas[$_] <= $intervalo)){
					$picos_intervalo{$masas[$_]} = $intensidades[$_];
					$num_peaks_intervalo++;
				}
			}
			my @masas_by_intensity;
			#ningun pico experimental en el intervalo
			if ($num_peaks_intervalo == 0){
				next;
			}
			#un solo pico experimental en el intervalo
			elsif ($num_peaks_intervalo == 1){
				foreach (keys %picos_intervalo ){
					push @masas_by_intensity, $_;
				}
			}
			#mas de 1 pico experimental en el intervalo
			else {
				foreach my $mass (sort {$picos_intervalo{$b} <=> $picos_intervalo{$a}} keys %picos_intervalo ){
					push @masas_by_intensity, $mass;
				}
			}

			#Matching entre masas teoricas y experimentales
			my $masas_testadas = 0;
			foreach my $masas_exp(@masas_by_intensity){
				$masas_testadas++;
				foreach (keys %fragmentos){
					if (abs($fragmentos{$_} - $masas_exp) <= $fragment_error){
						$correspondencias++;
						push @masas_exp_match_intervalo, $masas_exp;
						push @masas_teor_match_intervalo, $fragmentos{$_};
						push @tags_match_intervalo, $_;
					}
				}
				if($masas_testadas == $picos_x_sector){
					last;
				}
			}
		}
		# $num_fragments   => trials
		# $correspondencias =>  successes
		# $probablty  =>  numero picos_por_sector / 100
		my $probablty = $picos_x_sector / 100;
		my $binomial_acum_rango = binomial($num_fragments, $correspondencias, $probablty);
		my $score;
		if($binomial_acum_rango > 0){
			$score = -10 * (log($binomial_acum_rango) / log(10));
		}
		else{
			$score = 0;
		}
		#print "score, $score\n";
		$datos_num_picos{$picos_x_sector}{masas_exp_match_intervalo} = \@masas_exp_match_intervalo;
		$datos_num_picos{$picos_x_sector}{masas_teor_match_intervalo} = \@masas_teor_match_intervalo;
		$datos_num_picos{$picos_x_sector}{tags_match_intervalo} = \@tags_match_intervalo;
		$datos_num_picos{$picos_x_sector}{score} = $score;
	}
	return \%datos_num_picos;
}


1;      #787

__END__

=encoding utf8

=head1 NAME

=head1 SYNOPSIS

=head1 DESCRIPTION

=head1 LIMITATIONS

=head1 SEE ALSO

=cut
